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M.  V.  Mironenko 

report  on  the  project  R&D  8968-EN-01 
'DEVELOPMENT  OF  A  PROGRAM  TO  CALCULATE  LIQUID-LIQUID 
PHASE  EQUILIBRIA  IN  MULTICOMPONENT  SYSTEMS  CONSIST  OF 
ORGANIC  SUBSTANCES’. 

Abstract:  This  report  documents  a  Fortran  version  of  a  chemical  thermodynamic  model  for 
calculating  liquid-liquid  equilibria  in  mixtures  of  organic  compounds.  The  model  applies  the 
Gibbs  energy  minimization  method  for  phase  equilibria  computation  combined  with  the 
UN1FAC  routine  and  thermodynamic  database  for  calculating  activity  coefficients  of  organic 
substances  in  multicomponent  organic  liquids  on  the  basis  of  the  group  contribution  theory  The 
model  can  be  extended  without  modifications  of  the  Gibbs  energy  minimization  program  to  take 
into  account  also  chemical  interactions  between  organic  components. 

1 .  Statements  of  the  task. 

A  task  to  calculate  a  phase  composition  of  multicomponent  organic  liquid  mixtures  is  a  part  of 
more  general  task  to  calculate  chemical  equilibrium  composition  of  multicomponent  system. 
Because  of  this,  in  spite  of  the  components  at  specified  conditions  are  considered  as  inert  ones 
(in  other  words  they  do  not  interact  each  with  other),  to  keep  generality  we  will  use  the  terms  of 
chemical  equilibrium. 

We  have  a  system  consists  of  JV  organic  compounds  (components)  at  specified  temperature  and 
pressure.  The  total  amount  of  every  rth  compound  is  equal  to  Bj.  It  takes  to  determine  a  number  P 
of  the  phases  in  the  system,  concentrations  Xy  of  any  rth  component  in  every  yth  phase,  and 
masses  of  the  phases  Xj  (j=l ,P)  in  the  system.  It  is  assumed  that  any  component  is  present  in 
every  coexisting  phase. 

2.  _ Thermodynamic  conditions  of  multicomponent  equilibria. 

The  necessary  and  sufficient  condition  of  multicomponent  heterogeneous  equilibrium  is  that 
the  Gibbs  energy  of  the  system  G  is  minimum  under  the  balance  constrains.  The  Gibbs  energy 
function  of  the  system  is  as  follow: 

~  &  ~  "*■  Iri(^y))  (l) 
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where  is  activity,  is  the  amount,  and  ^  is  the  standard  chemical  potential  of  /th 
component  in  yth  phase.  If  we  consider  the  components  as  inert  ones  and  that  the  components 
have  the  same  standard  potential  in  any  phase,  the  values  of  zero  can  be  assigned  to  these 
parameters. 

nij  nij 

Activity  of  the  components  can  be  expressed  as  a <y  =  ,  where  x(j  =— f-  =  y —  IS 

j  Z  n,j 

/=  i 

a  mole  fraction  of  the  component  /  in  the  phase  j,  y}  is  activity  coefficient. 

The  consequence  of  the  condition  of  a  thermodynamic  equilibrium  is  that  the  chemical 
potential  for  each  /th  component  must  be  the  same  in  any /th  phase.  If  every  component  has  the 
same  standard  potential  ju  in  any  phase  it  leads  to  equality  of  activities: 

xi\Y<l  =  xi’2/i2  =  —=  *iPYiP  (2) 

Mass  balance  constrains  may  be  written  as  a  system  of  liner  equations; 

Z  n..  -Vjj  =  Q,i  =  \N  (3) 

J- 1 

where  Vy  is  a  stoichiometric  coefficient,  that  shows  a  number  of  independent  components  J 
in  the  component  i.  If  we  choose  as  independent  components  real  components  of  the  system  then 

yij  =  • 

It  can  be  shown  by  liner  transformations  that  P  s  N  (The  Gibbs  phase  rule) 

Summary  of  the  algorithm  of  the  Gibbs  energy  minimization. 

The  applied  algorithm  of  search  of  the  global  minimum  of  the  function  (1)  under  constrains 
(3)  combines  methods  of  linear  and  nonlinear  programming  It  is  very  close  to  the  algorithm 
applied  by  de  Capitani  and  Brown,  1987 
The  initial  phase  assemblage. 

We  start  with  a  feasible  assemblage  of  N  linearly  independent  pure  phases  consisted  of  only 
one  component  each.  The  amount  of  each  /th  phase  is  equal  to  B,  ,  Thus  the  mass  balance 
equations  are  satisfied. 

Step  1.  The  nonlinear  programming 

At  the  first  step  we  suppose  that  every ytb  potential  phase  has  the  following  composition; 
*,,=0.95  if  /=/  and  x,;=(l-0  95)/(7V-/y  if  /#,  The  activity  coefficients  yjof  the  components  in 
every  J  th  phase  can  be  calculated  with  the  UNIFAC  routine  for  these  specified  compositions  of 
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the  phases,  The  meaning  of  this  key  action  is  to  assign  different  initial  approximations  of  activity 
coefficients  to  the  components  in  various  potential  phases.  This  is  the  way  to  start  from 
extremely  different  meanings  of  concentrations  of  the  components  in  different  phases  that  have 
to  give  the  same  values  of  activities  to  be  produced  at  every  iteration  (see  equation  (2)). 

A  logarithm  of  the  equilibrium  constant  of  reaction  of  a  given  component  formed  from 
independent  components  of  the  system  is  equail  to  —AG  of  the  reaction, 

N 

In au  =  ln(x;J  Yu)  =  -AG  =  Vjj  ~ 

j 

It  means  that  for  every  potential  phase  J  the  all  mole  fractions  xu  arc  calculated  as 


N 


xij  =exPj  -<Mi  -ZMjVy  +Yu ) 


N 


XjJ 


The  exponent  is  always  >  0,  thus  we  can  divide  all  xf  by  their  sum  Xj  =  £  Xjj  :  =- 

/=l  Xj 


Now  the  sum  ofx,-  is  equal  to  1 

The  stoichiometric  vector  of  the  phase  J  can  be  written  as  vu  =  vu  -  xu .  The  value  of 
N 

-In  Xj  +  Zfifvjj  is  assigned  to  jjj 
/=  1 

At  the  next  step  these  phases  will  be  considered  as  having  a  fixed  composition. 

Step  2.  Linear  programming 

The  problem  to  find  the  equilibrium  assemblage  from  a  set  of  fixed  composition  phases  can  be 
formulated  as: 

M 

minimize  G  =  £  rtj >  0  under  constraints  (3)  ,  M=2N  is  a  number  of  the  independent 

components  plus  a  number  of  potential  phases. 

This  can  be  solved  with  a  simplex  method.  Simplex  is  a  classic  finite  iteration  method  of  linear 
programming  (Korn  and  Korn,  1963).  The  following  procedure  assumes  that  at  the  beginning  of 
step  2  the  first  N  columns  in  matrix  Jv|  define  a  feasible  assemblage  of  linearly  independent 

phases  with  and  masses  B\ ,  These  conditions  are  satisfied  by  the  initial  assemblage  as  defined 
earlier  and  remain  so  throughout  all  iterations.  The  matrix  ||vj|  is  as  follow: 


1  0  0...0  V ifN+i  Vi?n+2  . r.Vi^N  Bj 
0  1  0...0  V2,N+I  V2,1SH-2  . .  -  V2,2N  B2 

0  0  0...1  Vn,nh  VNfN,2  ...VN,2N  Bn 


3 


'30/03  '01  12:51  FAX  9382054 


GEOHI  RAN 


@05 


Each  column./  =  1 ,2. . .  N  represents  a  “currently  stable  phase.”  Each  column  A=N+1 ,  N+2, 

.  --2N  can  be  interpreted  in  terms  of  a  reaction  relating  the  phase  k  to  the  "currently  stable 
phases".  Gibbs  energy  of  the  reaction  is 
N 

AGk^Mk  ~  £  MjVjk 
/-I  1 

Is.  Each  reaction  is  tested  for  G(k)  <  0.  If  one  is  found  then  that  phase  k  is  more  stable  than  a 
linear  combination  of  the  “currently  stable"  phases.  At  this  point  the  following  procedure  is 
started: 

2s.  Locate  the  "currently  stable  phase"  to  be  replaced  by  the  newly  found  phase  k.  The 

smallest  positive  B/vJk is  the  criterion.  Phase  j  will  be  replaced  by  phase  k.  Thermodynamic 

meaning  of  this  is  the  reaction  of  new  phase  formation  can  go  until  one  of  the  reactants  is  used 
up 

3  s.  Replace  phase  j  by  phase  k. 

4s.  Reduce  [|vJJ  This  normalizes  the  reactions  to  refer  to  the  phases  now  occupying  the  basis 
(the  first  N  positions  of  the  matrix) 

Return  to  step  1  s.  If  none  is  found  terminate  the  linear  programming. 

Return  to  step  1  (nonlinear  programming),  Repeat  in  series  steps  1  and  2  until  for  every 
Xj  in  every  potential  phase  the  condition  |l  -  x(/*V  x(/)|  <  £  fulfils,  (j)  is  a  number  of  iteration 
and  e  is  a  needed  accuracy. 

As  a  result  of  the  iterative  procedure  the  basis  of  the  matrix  will  be  occupied  by  stable 
phases.  Some  of  this  (or  all,  if  the  solution  is  a  homogeneous  system)  will  have  the  same 

composition.  Such  phases  are  considered  as  parts  of  one  phase,  and  their  amounts  Bj  should  be 
summarized. 

Activity  coefficients  calculation. 

To  calculate  activity  coefficients  of  compounds  in  liqiud  multicomponent  phases  as 
function  of  composition  we  used  UNEFAC  prediction  computer  program  and  database 
(Fredenslund  et  al.,  1974).  The  group  contribution  method  is  based  on  the  following 
fundamental  assumptions.  1)  The  logarithm  of  activity  coefficient  is  a  sum  of  two 
contributions:  a  combinatorial  part,  due  to  differences  in  size  and  shape  of  the  molecules,  and 
the  residual  part,  due  to  energy  interactions.  2)  The  residual  part,  the  contribution  from  group 
intei  action,  is  assumed  to  be  the  sum  of  the  individual  contributions  of  each  solute  group  in  the 
solution  less  the  sum  of  the  individual  contributions  in  the  pure-component  environment,  3) 

The  individual  group  contributions  in  any  environment  containing  groups  of  kinds  J,  2  A' are 
asuumed  to  be  only  a  function  of  group  concentrations  and  temperature  (Van  Ness  ct  al, ,  1973, 
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Reid  ct  al.,  1987).  The  database  on  group  interaction  parameters  had  been  developed  as  a  result 
of  processing  of  large  amount  of  experimental  data. 

The  program  is  written  in  FORTRAN-IV  and  consists  of  the  following  modules. 

1.  Main  program  EQSIMFR  minimises  Gibbs  energy  of  the  system  and  outputs  the  result 
of  calculations  It  calls  UNIFAC  subroutine  first  time  to  read  input  data  and  to  calculate 
parameters  for  specified  system  and  N  (number  of  components)  times  at  every  iteration  to 
calculate  activity  coefficients  of  the  components  in  potential  phases  at  their  current 
concentrations.  The  program  also  calls  SIMPLEX  subroutine  to  search  the  assemblage  of  fixed 
composition  phases  at  every  iteration. 

We  used  UNIFAC  program  and  database  of  Fredenslund  ct  al.  (1992)  without  significant 
changes.  The  only  corrections  were  to  convert  the  program  to  a  subroutine  and  to  separate  the 
part  of  the  routine,  wh  ich  searches  in  the  database  and  calculates  parameters  for  specified 
system  to  execute  it  one  time. 


We  kept  input  without  any  changes.  The  information  about  the  phase  composition  in 
Fredenslund  et  al.  program  is  used  in  the  model  as  a  bulk  composition  of  the  system 

An  example  of  the  inmitffile  “input.txf’Y 


(1)  PROP.  GLYCOL 

(2)  ETHANOL 

(3)  WATER  1 

(4)  BENZENE  1 


4  Name  of  the  componenl(20  positions)  and 
3  a  number  N  of  different  secondary  groups  in  the  molecule 


1  1  3 

1  1  2 
1  17 

6  9 

302.50 
0.005 


2  1  3  2  14  Number  of  limes  secondary  group  I  appears  in  molecule  J, 

2  1  14  Identification  number  of  group  /  in  molecule  J,  secondary 

group  number 

temperature 

0.06  0.935  0. 1  numbers  of  moles  (or  mole  fractions)  of  the  components 

in  mixture 


An  example  of  the  output  (file  “result”! 


TEMPERATURE  302.50  K 
PHASE  1 
0.1694  MOLES 

n  COMPONENT  MOLE  MOL.FRACTION  ACTIVITY  ACT.COEFF. 
1  (1)  PROP.  GLYCOL  0.50006D-02  0.29514D-01  0.33594D-01  1.138 
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2  (2)  ETHANOL  0.59736D-01  0,352571>K)0  0.38978D+00  1.106 

3  (3)  WATER  0.473 11D-02  0.27923D-01  0,99966D+00  35.800 

4  (4)  BENZENE  0.99964D-01  0.590000+00  0.59293D+00  1.005 

PHASE  2 

0.9306  MOLES 

n  COMPONENT  MOLE  MOL.FRACTION  ACTIVITY  ACT.COEFF. 

1  (1)  PROP.  GLYCOL  0.31332D-06  0.33670D-06  0.33594D-01  99773.219 

2  (2)  ETHANOL  0.26361D-03  0.28328D-03  0.38978D+00  1375.930 

3  (3)  WATER  0.93027D+00  0.99968D+00  0.99968D+00  1.000 

4  (4)  BEN2IENE  0.3G295D-04  0.39003D-04  0,592930+00  15202.278 

BALANCE  (MOLES) 

COMPONENT  PHASE  1  PHASE  2  BULK  CALC.  BULK  INPUT 

(1)  PROP.  GLYCOL  0.0049997  0.0000003  0.0050000  0.0050000 

(2)  ETHANOL  0,0597364  0.0002636  0.0600000  0.0600000 

(3)  WATER  0,0047313  0.9302687  0.9350000  0,9350000 

(4)  BENZENE  0.0999637  0.0000363  0.1000000  0,1000000 

Concluding  remarks 

It  was  noted  that  the  model  could  also  take  into  account  chemical  interactions  between 
substances.  For  this  goal  it  takes  I)  input  standard  chemical  potentials  fJ°  of  the  substances  at 
specified  temperature  and  2)  the  matrix  |jv]| contains  real  stoichiometric  (in  chemical  element  or 
in  group  terms)  compositions  of  the  substances. 
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